setup;
global h
global k

syms sim_x sim_y sim_u sim_f;
% sim_u = 7*sim_x^7*sim_y^7+2*sim_x^3*sim_y^3+3*sim_x^9*sim_y^2+8*sim_x^6*sim_y^2;
% sim_f = diff(sim_u,sim_x,2) + diff(sim_u,sim_y,2);

% 
% 
% x = 2 : h : 14;
% y = 3 : h : 16;
% 
% m = length(x);
% n = length(y);
% % 
% u = reshape(subs(subs(sim_u, x), y), m,n);
% f = reshape(subs(subs(sim_f, x), y), m,n);
% 

m = 75;
n = 36;
h = 1/(m-1);
x = 0:h:1;
y = 0:h:1;
u = zeros(m,n);
f = u;
W = u;
epsilon = 1;

for i = 1:m
    for j = 1:n
        u(i,j) = 3*x(i)*y(j)^2+ x(i)^3 + 2*y(j)^3 + y(j)^4; 
        %f(i,j) = 6*x(i) + 6*x(i) +12*y(j) + 12*y(j)^2; 
        f(i,j) = 6*x(i)*(1+sin(x(i))+cos(y(j)))+3*(x(i)^2+y(j)^2)*cos(x(i))-2*y(j)*(3*x(i)+y(j)*(2*y(j)+3))*sin(y(j))+2*y(j)*(4*y(j)+3)*(sin(x(i))+cos(y(j))+1)+2*(3*x(i)+y(j)*(2*y(j)+3))*(sin(x(i))+cos(y(j))+1);
        W(i,j) = 1 + epsilon .* (sin(x(i)) + cos(y(j)));
    end
end

k = floor(log2(min(m,n)))-1;

% u = zeros(m,n);
% f = rand(m,n);

g=zeros(m,n);
g(1,:) = u(1,:);
g(:,1) = u(:,1);
g(m,:) = u(m,:);
g(:,n) = u(:,n);

tic
v = solverPoi(crop(f), g, W);
toc
figure(3)
surf(abs(crop(u) - v));

u_c = crop(u);
norm(u_c(:) - v(:))
